Thermodynamics of two lattice ice models in three dimensions 
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In a recent paper we introduced two Potts- like models in three dimensions, which share the follow- 
ing properties: (A) One of the ice rules is always fulfilled (in particular also at infinite temperature, 
/S = 0). (B) Both ice rules hold for groundstate configurations. This allowed for an efficient calcula- 
tion of the residual entropy of ice I (ordinary ice) by means of multicanonical simulations. Here we 
present the thermodynamics of these models. Despite their similarities with Potts models, no sign 
of a disorder-order phase transition is found. 

PACS numbers: PACS: 05.50.+q, 11.15.Ha, 12.38.Gc, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 



By experimental discovery [l[ it was found that ice I 
(ordinary ice) has in the zero temperature limit a resid- 
ual entropy S/N = k ln(Wi) > where N is the number 
of molecules and W% the number of configurations per 
molecule. Subsequently Linus Pauling Q based the esti- 



mate W l 



Pauling 



3/2 on the ice rules: 



1. There is one hydrogen atom on each bond (then 
called hydrogen bond). 

2. There are two hydrogen atoms near each oxy- 
gen atom (these three atoms constitute a water 
molecule). 

Pauling's combinatorial estimate turned out to be in ex- 
cellent agreement with subsequent refined experimental 
measurements Q. This may be a reason why it took 
25 years until Onsager and Dupuis Q pointed out that 
W\ = 1.5 is only a lower bound. Subsequently Nagle @ 
used a series expansion method to derive the estimate 

^Naglc = L5f j 685 ( 15 ) 

In Q we introduced two simple models with nearest 
neighbor interactions on 3D hexagonal lattices, which al- 
low one to calculate the residual entropy of ice I by means 
of multicanonical (MUCA) @, 1, ij simulations. The 
hexagonal lattice structure is depicted in Fig. [T] [l(| • In 
the first model, called 6-state (6-s) H2O molecule model, 
ice rule (2) is always enforced and we allow for six dis- 
tinct orientations of each H2O molecule. Its energy is 
defined by 
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FIG. 1: Lattice structure of one layer of ice I (reproduced 
from Ref. [7]). The up (u) sites are at z = l/v24 and the 
down (d) sites at z = — l/\f2A. Three of its four pointers to 
nearest neighbor sites are shown. 



Here, the sum is over all bonds b of the lattice (sj and 
si indicate the dependence on the states of the two H2O 
molecules, which are connected by the bond) and 



h(h s 1 <* 2 \ - / 1 for a n y dro S en bond : 
W>*b> 8 b)- to otherwise. 



(2) 



In the second model, called 2-state (2-s) H-bond model, 
ice rule (1) is always enforced and we allow for two posi- 
tions of each hydrogen nucleus on a bond. The energy is 
defined by 



E 



E 



f(s,blb 2 s ,b 3 s ,bt) 



(3) 



where the sum is over all sites (oxygen atoms) of the 
lattice and / is given by 



* Corresponding author. 



f(8,b\,blb 3 a ,bi) 



(4) 
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{2 for two hydrogen nuclei close to s , 
1 for one or three hydrogen nuclei close to s , 
for zero or four hydrogen nuclei close to s . 

The groundstates of either model fulfill both ice rules. 

In this paper we use units with k = 1 for the Boltz- 
mann constant, i.e., (3 — 1/T. At (3 = the number of 
configurations is 6^ for the 6-s model and 2 2N for the 
2-s model. This sets the normalization, which can then 
be connected by a MUCA simulation of the type [8j to 
(3 values large enough so that groundstates get sampled. 
In reasonably good agreement with Nagle the estimate 
w mvga = L50738 ( 16 ) was obtained in @. In Ref. [U 

these calculation were extended to partially ordered ice 
for which corrections to groundstate entropy estimates 
by Pauling's method were previously not available. 

A considerable literature [H [H, [jj, [H [H, [TJ H 03 
exists on lattice ice models. Most of these papers deal 
with 2D square ice. An extension to 3D is considered 
in [IH, [3, [HI] ■ All these models have in common that 
they enforce both ice rules generically and not just for 
the groundstates. So, they are non-trivial at all coupling 
constant values, while it is precisely the triviality of our 
ice models at (3 = 0, which allows one to set the normal- 
ization for the entropy and free energy, and they share 
with certain spin models [2(| that the residual entropy 
of their groundstates violates the third law of thermody- 
namics. 

Superficially our models are similar to g-state Potts 
models [2l| with q = 6 and the Ising case q = 2, which 
have first (q = 6) and second (q = 2) order phase transi- 
tions in 2D as well as in 3D. In contrast to that, we pro- 
vide numerical evidence in this paper that our ice models 
do not undergo a disorder-order phase transition at any 
finite value of (3. Our results are presented in section HT1 
Summary and conclusions follow in section ITTT1 
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FIG. 2: Energy per site for the 6-s and 2-s models. 



tice sizes used are compiled in table |TJ The lattice con- 
tains then N — n x n y n z sites, where n x , n y , and n z are 
the number of sites along the x, y, and z axes, respec- 
tively. Periodic BCs restrict the allowed values of n x , 
n yi and n z to n x = 1, 2, 3, . . ., n y = 4, 8, 12, . . ., and 
n z = 2, 4, 6, ... . Otherwise the geometry does not close 
properly. 

As proposed in [23| we used a Wang-Landau [24| recur- 
sion for determining the MUCA weights and performed 
subsequent MUCA data production with fixed weights. 
With one exception we used 32 x (20 x 10 6 ) sweeps per 
lattice for data production. For the largest lattice of the 
6-s model we produced a ten times larger statistics. Ta- 
ble U listed for each lattice size and model the number of 
cycling events from the average disordered energy Eq at 
(3 = to the groundstate energy E g and back, 



En 



E, 



II. SIMULATION RESULTS 



TABLE I: Overview of our multicanonical simulations. Here, 
, i^y , Hz are the number of lattice sites along the x, y, z 
axes, and N — n x n y n z . 



9 • 



(5) 



III 


n y 


n z 


N 


cycles (6-s) 


cycles (2-s) 


4 


8 


4 


128 


37828 


141825 
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12 
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288 


9455 


33205 
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12 
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360 


4891 


21621 
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12 
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576 


653 


11479 
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16 


8 


896 


412 


6452 


8 


20 


10 


1600 


215 


1587 


10 


24 


12 


2880 


1133 


506 



Using periodic boundary conditions (BCs), our sim- 
ulations are based on a lattice construction 22] simi- 
lar to that set up for Potts models in 9]. The lat- 



as recorded during the production part of the run. From 
the energy functions (fTJ) and ([3]) one finds Eq — —N for 
the 6-s model (there are two hydrogen atoms per oxygen 
and the probability to form a hydrogen bond is 1/2), 
Eq = — 1.25 N for the 2-s model (at one site there are 
16 arrangements of hydrogen atoms with average energy 
contribution -[2x0 + 8xl + 6x 2]/16 = -1.25), and 
Eg = —2N for both models. In the following we restrict 
the (3 range of our figures to < (3 < 5, which is large 
enough to sample groundstates in sufficient numbers so 
that extrapolations down to temperature T ~ become 
controlled. 

In Fig. [2] we show the average energy per site, E/N, 
from the MUCA simulations of our two models as ob- 
tained by the reweighting procedure [9J (note that we 
use E for the energy of configuration as well as for av- 
erage values over configuration energies and assume the 
reader knows to distinguish them). Obviously there are 
almost no finite size effects, because the curves from all 
lattice sizes fall within small statistical errors, which are 
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FIG. 3: Specific heat per site for the 6-s and 2-s models. FIG. 4: Specific heat per site for the 2D Ising model on JV = 

L 2 lattices. 



not visible on the scale of this figure, on top of one an- 
other. 



TABLE II: Some specific heat data C/N with error bars (in 
parenthesis) for the N = 2880 lattice. 



(3 


6-s model 


2-s model 


0.5 


0.1175119 (77) 


0.093780 (17) 


1.5 


0.673681 (71) 


0.48331 (11) 


2.5 


0.87873 (19) 


0.59913 (22) 


3.5 


0.69110 (35) 


0.46235 (42) 


4.5 


0.43066 (45) 


0.28637 (61) 



The specific heat per site, C/N, is calculated via the 
fluctuation-dissipation theorem, 



(6) 



and plotted in Fig. [3l Finite size corrections are now 
visible for the smallest, N = 128, lattice. For the other 
lattices the curves fall within error bars on top of one 
another. Error bars were calculated with respect to 32 
jackknife bins and are at some j3 values included for our 
largest, N = 2880, lattice. Some data for these points are 
given in table HH Note that the N = 2880 data for the 
6-s model rely on a ten times large statistics than those 
for the 2-s model, while the error bars are only slightly 
smaller. As noticed before 0], the simulations of the 2-s 
model are more efficient for determining the groundstate 
entropy than simulations of the 6-s model. Fluctuations 
increase with lattice size, so that it is more difficult to 
obtain accurate results on large than on small lattices. 

We want to contrast Fig. [3] with specific heat results for 
the 6-state and 2-state Potts models on L D lattices. Im- 
mediately, one notices that it is not entirely clear whether 
this comparison should be done in 2D or 3D. While the 
space dimension in which our ice models are embedded 
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FIG. 5: Specific heat per site for the 3D 6-state Potts model 
on N = L 3 lattices. 



is clearly 3D, each site is connected through links with 
four neighboring sites, which is the Potts model situation 
in 2D. The 2D and 3D Ising models are well known for 
their second order phase transitions. The specific heat 
is logarithmically divergent in 2D [2|| and has a critical 
exponent awO.l in 3D (see [26| for a review). The 2D 
and 3D 6-state Potts models have first order transitions 
with a larger latent heat per spin in 3D than in 2D (in 
the normalization of [j| AE/N = 0.40292828 in 2D (28| 
and AE/N = 2.36442 ± 0.00017 in 3D [H). 

For second order transitions the maximum of the spe- 
cific heat diverges ~ ln(L) for a logarithmic divergence 
(q = 0) and ~ L"'" for a > 0, where v is the criti- 
cal exponent of the correlation length. In case of first 
order phase transitions the peak in the specific heat di- 
verges ~ L D , where the proportionality factor is [27| 
(f3t) 2 (AE / N) 2 with (3 t the inverse transition tempera- 
ture and AE/N the latent heat per spin. 
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N=128 
N=288 
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N=576 
N=896 
N=1600 
N=2880 



2-s model. Relative statistical errors are smaller than 
those in Fig. [3] for the specific heat. In the (3 — » oo 
limit our data improve slightly on the results reported 
in Ref. Q, because we have with N = 2880 one larger 
lattice added. Consistent fits to the previously discussed 
form Wi(x) = W^ VGA + a lX e , x=l/N combine to 



MUCA 



1.50721 (13) and = 0.901(16) 



(7) 



The error bars in parenthesis are purely statistical and 
do not reflect eventual, additional systematic errors due 
to higher order finite size corrections. 
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FIG. 6: Free energy 



III. SUMMARY AND CONCLUSIONS 
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FIG. 7: Entropy . 



In Figs. |4] and [5] we plot the specific heat on various 
lattices for the two extremes, the weak logarithmic diver- 
gence for the 2D Ising model and the strong divergence 
for the 3D 6-state Potts model. For the 2D Isin g m odel 
the analytical solutions of Ferdinand and Fisher [3CJ] are 
plotted, while the plots for the 3D 6-state Potts model 
rely on recent numerical results [29j . It is clear that even 
the case of a weak logarithmic divergence is markedly 
distinct from the behaviors in Fig. [3j where no finite size 
effects are observed within the rather accurate statistical 
errors. This distinction becomes all too obvious when 
the comparison is made with the strong first order phase 
transition of the 3D 6-state Potts model. 

To complete the picture of our two ice models we plot 
in Figs. [HI and [7] their free energy and entropy densities as 
obtained from our simulations, using as input the known 
normalizations at f3 — 0. In the cases at hand these 
are S /N = ln(6) for the 6-s and S /N = ln(4) for the 



The unusual properties of water and ice owe their exis- 
tence to a combination of strong directional polar inter- 
actions and a network of specifically arranged hydrogen 
bonds [3l|, [H, [33| . The groundstate structure of such a 
network can be described by simple lattice models, which 
are defined by the energy functions (JTJ) and ([3]). 

In the present paper we have presented finite size scal- 
ing evidence that there is no phase transition between 
(3 = and the groundstate region of these models. This 
lack of a transition makes reliable estimates of the com- 
binatorial groundstate entropy of ice I particularly easy. 

Tentatively, we like to see a reason for the marked 
difference to q = 6 and q = 2 Potts models in the 
large groundstate entropy, S/N = ln(Wi), of our ice 
models, which violates the third law of thermodynamics, 
while the groundstate entropy of Potts models, S/N = 
\n(q)/N, approaches zero in the N — > oo limit. This 
is not an entirely convincing argument as the effective 
number of states W\ per spin is still about 2 (i.e., larger 
than 1.5) for the 3D 6-state Potts model at the ordered 
endpoint of the transition [2i| . 

Note that we did not investigate bond statistics in the 
groundstate ensemble, which one may expect to exhibit 
critical correlations. 
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